History-independence of steady state in simultaneous two-phase flow through porous 

media 



m 

o 

(N 

u . 

< 



q 

O 



m 

o 
en 



% 



Marion Erpelding,^ Santanu Sinha,^ Ken Tore Tallakstad,'^ 
Alex Hansen,^ Eirik Grudc Flckk0y/ and Knut J0rgen Mal0y^ 

^Department of Physies, University of Oslo, PB 1048 Blmdern, N-0316 Oslo, Norway 
'^Department of Physics, Norwegian University of Science and Technology, N-7491 Trondheim, Norway 

(Dated: 18 avril 2013) 

It is well known that the transient behavior during drainage or imbibition in multiphase flow in 
porous media strongly depends on the history and initial condition of the system. However, when 
the steady-state regime is reached and both drainage and imbibition take place at the pore level, 
the influence of the evolution history and initial preparation is an open question. Here, we present 
an extensive experimental and numerical work investigating the history dependence of simultaneous 
steady-state two-phase flow through porous media. Our experimental system consists of a Hele- 
Shaw cell flUed with glass beads which we model numerically by a network of disordered pores 
transporting two immiscible fluids. From the measurements of global pressure evolution, histogram 
of saturation and cluster-size distributions, we find that when both phases are flowing through the 
porous medium, the steady state does not depend on the initial preparation of the system or on the 
way it has been reached. 



I. INTRODUCTION 

Understanding the physical mechanisms underlying 
multiphase flow in porous media is crucial to a wide va- 
riety of industrial and environmental problems, such as 
oil recovery, CO2 transport and storage or ground water 
management. At the macroscopic scale, the flow of im- 
miscible fluids in a porous material is described by speci- 
fying the relations between global quantities such as flow 
rate, pressure gradient or fluid saturation. At the pore 
scale, this flow is governed by the competition between 
capillary, viscous and gravitational forces. Understanding 
the link between the two levels of description requires to 
relate the position and shape of the interface(s) between 
the two phases to the values of the macroscopic variables. 
From an experimental point of view, two-dimensional 
model porous media providing direct pore-scale visuali- 
zation of the flow structures are ideal tools to study mul- 
tiphase flow mechanisms and their relations with global 
quantities in controlled, laboratory-scale situations. Over 
the past decades, micro channel networks etched in trans- 
parent plates [1], or prepared using molding techniques 
[2] and porous Hele-Shaw cells consisting of a layer of 
beads between two parallel plates [3] have become clas- 
sical tools, to which improvements have been constantly 
proposed [4, 5]. Experimental observations have been ex- 
plained through extensive numerical simulations based 
on network models [6-10] and lattice Boltzmann methods 
[11-17], statistical models [18-20] and differential equa- 
tions [21]. Most of the research in this area has been 
focused upon transient phenomena, i.e. drainage or im- 
bibition - arising when one phase displaces the other in a 
porous medium. The relation between macroscopic flow 
parameters, fluid morphology and stability of the inter- 
face between the two phases has been thoroughly obser- 
ved [3, 4] and successfully modeled [18, 19]. 

In this article, we deal with steady-state flow, where 



many questions are yet to be answered. In the steady- 
state regime, two phases are injected simultaneously into 
the porous medium and one observes that one or both are 
fragmented and transported in the form of clusters of va- 
rious sizes, forming a complex flow pattern with multiple 
interfaces. After a characteristic time, the system reaches 
a steady-state in which the macroscopic flow variables 
fluctuate around constant values. The usual distinction 
between drainage and imbibition is irrelevant to describe 
steady-state two-phase flow, in which both processes oc- 
cur simultaneously. New approaches arc thus needed to 
understand this regime. In an effort to bring new insight, 
experimental and numerical studies have investigated the 
relations between macroscopic flow variables, and dif- 
ferent models have been proposed to relate them to pore- 
scale flow mechanisms : Payatakes and co-workers carried 
out detailed experimental, numerical and theoretical stu- 
dies of steady-state two-phase flow, with emphasis on the 
determination of relative permeabilities [1, 9, 10, 22-25]. 
The steady-state characteristics of macroscopic flow pro- 
perties have also been investigated numerically by Knud- 
sen et al. [26]. A power- law relation between pressure 
and steady-state flow rate has been observed experimen- 
tally by Tallakstad et al. in a 2D system [27, 28], and by 
Rassi et al. in a 3D system [29]. Very recently, the re- 
lation between the steady-state flow rate and pressure 
drop has been derived analytically for two-phase flow 
through single capillaries [30] and through porous me- 
dia [27, 28, 31] and also supported by extensive nume- 
rical simulation. Distributions of non- wetting or wetting 
clusters in the steady state have also been studied expe- 
rimentally by Tallakstad et al. [27, 28] and numerically 
by Ramstad and Hansen [32], and critical exponents were 
measured. It is worth noting though, that the compari- 
son of experimental and numerical results is not always 
straightforward, due to the different boundary conditions 
used in the two cases. Yet from the point of view of sta- 



tistical physics, the existence of a genuine steady state 
is very promising to build a thermodynamic-like theo- 
retical description of the system. In that context, it is 
crucial to determine whether the steady state is inde- 
pendent on the history of the process, or in other words, 
whether it is a real state in a thermodynamic sense [33]. 
It is well known that when drainage and imbibition oc- 
cur successively, the relative permeabilities become his- 
tory dependent and the pressure-saturation curves dis- 
play an hysteresis [34, 35], the underlying pore-scale me- 
chanisms of which are known. Besides in the well-known 
magnetic and elastic systems, such history dependence 
and hysteresis has also been observed in different flow 
processes like hydrodynamic heat flow [36] and particle 
flow through random media [37]. However, in the case 
of two-phase flow through porous media, it is not tri- 
vial to predict whether such an hysteresis will come into 
play in the steady-state situation when drainage and 
imbibition occur simultaneously. It was proposed that 
a thermodynamic-like description for simultaneous two- 
phase flow in porous media can be sketched [33] if the 
flow was history independent. 

Here, we present an extensive experimental and 
numerical study in order to investigate the history- 
independence of the steady state. Our experimental sys- 
tem consists of a porous Hele-shaw cell in which we simul- 
taneously inject air and a viscous water-glycerol solution. 
We then compare steady-states obtained for a given flow 
rate with different initial conditions. We model the sys- 
tem by a network of disordered pores transporting two 
immiscible fluids. From pressure measurements and ana- 
lysis of the statistical properties of the flow patterns, we 
observe no history-dependence for the steady state. 



II. EXPERIMENTAL SETUP 

We use a two-dimensional, transparent, porous Hele- 
shaw cell. This experimental setup, shown on Figures I 
and 2, has been described in details in [27] and we re- 
call its main features here. The porous medium consists 
of a random mono layer of glass beads, 1 mm in dia- 
meter, spread between the sticky sides of two sheets of 
contact paper. Its lateral boundaries are sealed with sili- 
con glue. Attached on top of this layer, a Plexiglas plate 
with etched flow channels allows injection and evacua- 
tion of fluids into and from the porous matrix. A pres- 
sure cushion (see [27] for details) placed below the porous 
medium, and a thick glass plate on top prevent the sys- 
tem from bending when the pressure increases as fluids 
are injected. Clamps placed all around the setup main- 
tain all the layers together. This way, we obtain a po- 
rous medium of constant thickness a — 1 mm, length 
L = 85 cm, width W = 42 cm in which the beads remain 
immobile. The porosity 4> and absolute permeability kq 
of the medium are found experimentally to be </> = 0.63 
and Ko = 1.95- 10^ cm^ [27]. 

We use the same fluid pair as Tallakstad et al. [27], 
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Figure 1. Sketch of the experimental setup : the 2D po- 
rous matrix consists of a disordered mono layer of glass beads 
spread between two sheets of contact paper. The boundaries 
are sealed with silicon glue. The upper part of the system 
consists of a Plexiglas plate with drilled inlet and outlet flow 
channels. A pressure cushion and a thick glass plate placed 
below and above the porous matrix ensure the overall rigidity 
of the system and maintain its thickness constant. Clamps 
maintain all the layers together. A lightbox illuminates the 
system from below and a digital camera is placed above to 
record images of the flow structure. 



namely air as the non- wetting phase and a viscous water- 
glycerol solution (15-85% in mass) as the wetting phase. 
The latter is dyed in black with 0.1% Ncgrosine to ob- 
tain a good contrast on the experimental images (see 
Figure 3). As pointed out by Tallakstad et ai, the use 
of air as one of the phases introduces more complexity 
in the two-phase flow problem for two reasons : flrst, it 
yields a high viscosity contrast with the glycerol solu- 
tion, and second, its compressibility gives rise to rapid 
bursts or avalanches [27]. However, from an experimen- 
tal point of view, it has the huge advantage of allowing 
us to reuse the same porous model for all experiments. 
Indeed, it can easily be flushed out, making it possible 
to obtain reproducible initial conditions. Therefore, we 
find it fully convenient for the present study. As shown 
on Figure 2, the two phases are injected simultaneously, 
using the same syringe pump, from 15 syringes (7 of air 
and 8 of water-glycerol solution) , each connected to one 
of the 15 inlet nodes of the porous model. In the follo- 
wing, we note Q the total flow rate, while Qv, — {8/l5)Q 
and Qnw = (J/15)Q denote the wetting and non-wetting 
flow rates, respectively. To account for small temperature 
variations due to the heat released by the lightbox (see 
Figure 1), we monitor the temperature of the wetting 
phase at the outlet of the model. The viscosity /!„, of the 
wetting phase is deduced accordingly using the empirical 



Syringe pump 



Imposed flow rate 



\/ 

Pressure 
sensor (P.S) 




' Thermometer 



Figure 2. Sketch of the experimental setup with the injection 
system. The two phases are contained in 15 syringes, each 
connected to one of the 15 inlet nodes of the porous model (7 
syringes of air represented in white and 8 syringes of water- 
glycerol solution in black). The same syringe pump is used to 
inject both phases simultaneously. The dotted lines give the 
dimensions of the area studied by image analysis (note that 
proportions are not respected). 



formula given in [38]. In the series of experiments pre- 
sented here, the measured temperatures are in the range 
24.4 - 29.3° Celsius, giving 0.083 > ^i.,^ > 0.062 Pa.s. 
The viscosity of air being fin^ w 1.9 x 10~'^ Pa.s, the 



viscosity ratio M = ^„ 
experiments. 



j/Hvj of the order of 10 in aU 



With this setup, the total flow rate Q and the fractio- 
nal flow Fw = Qw/Q are controlled flow variables. The 
volumes of wetting and non- wetting fluids present in the 
porous matrix, V^, and Vnwj are free to vary with time. 
Thus, the saturations S^ = V^jV and ^nw = Kiw/^j 
where V is the total pore volume, are free to fluctuate. 
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Figure 3. (a) Example of steady-state image (2400 x 3800 
pixels) : the flow pattern is made of air clusters (gray) of 
various sizes surrounded by the viscous liquid (black) - (b) 
Zoom : the high image resolution makes it possible to distin- 
guish glass beads (bright gray), air clusters (gray) and viscous 
liquid (black). These three phases yield three "peaks" on the 
gray scale image histograms, as illustrated by (c) . The height 
of the peaks contains information about the saturation of the 
system. 



The flow is characterized by the capillary number : 



Ca: 






(1) 



where fi^ is the wetting phase viscosity, Q^ is the 
total wetting fluid flow rate, 7 « 6.4 • 10~^ N/m is 
the interfacial tension between the two phases [27] and 
A = Wa(j> is the cross-section of the porous matrix. In 
the present experiments, we have explored the range 
3.33 10-s ^ Ca ^ 1.13 10"''. The highest experimental 
value of Ca is set by the maximum pressure that the 
porous model can hold. However, as we will see in Section 
IV, we have also explored higher values of Ca in numeri- 
cal simulations, namely 1.92 x 10~^ ^ Ca ^ 7.0 x 10^^. 

Our analysis of steady-state flow and its history depen- 
dence relies on two kinds of information : measurements 
of the pressure inside the model, and pictures of the flow 
pattern. We measure the pressure P{t) in the wetting 
phase as a function of time t using flow-through pressure 
sensors (Sensor Technics 26PC0100G6G) placed at three 
different points of the model as indicated on Figure 2. 
The porous model is lit from below using a lightbox, and 
images of the flow structure are recorded regularly using 
a Nikon D200 digital reflex camera giving 2592 x 3872 
pixels images with a spatial resolution of 8 x 8 pixels per 
mm^. Before further processing, images are cropped to 



remove boundaries, leaving us with 2400 x 3800 pixels on 
the final images. Figure 3(a) shows an example of steady- 
state image. The area of the imaged zone (see Figure 2) is 
large enough to contain many air clusters of various sizes. 
As illustrated by Figure 3(b), the high image resolution 
makes it possible to distinguish glass beads, air clusters, 
and viscous liquid. Each of these phases gives rise to a 
peak on the gray-scale image histogram (see Figure 3(c)). 
The heights of these peaks contain information about the 
proportion of wetting and non- wetting fluids in the sys- 
tem, thus about the saturations S'w and S„^. 



III. EXPERIMENTAL RESULTS AND 
DISCUSSION 

A. Principle of the experiments 

We investigate the history-dependence of steady-state 
flow by comparing steady states obtained at the same 
flow rate but with different initial conditions. For this, 
we have performed experiments in which the flow rate is 
modified twice, as illustrated by Figure 4(a). The porous 
model is initially filled with the wetting phase only. Then, 
both phases are injected simultaneously at a fixed flow 
rate Qi- Once the system has reached a steady state (ssi), 
we abruptly change the flow rate to a different value Q2 
and wait until a new steady state is established (SS2). 
Finally, we change the flow rate back to its initial value 
Qi and let the system evolve towards a third steady state 
(sss). We have used different couples {Qi,Q2}, listed in 
Table I, to obtain different magnitudes and signs of AQ = 
Qi — Qi- Wc compare ssi and SS3 using three criteria : 
the average pressure drop across the system, the fluid 
saturations, and the size distribution of the air clusters. 
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Table I. Experimental parameters and corresponding mea- 
surements. Qi and Q2 refer to the total imposed flow rates 
{i.e., for 15 syringes). Xiniet and Xoutiet are computed from 
the measured pressure drops according to Eq. 3. Note that 
experiments Exp.l and Exp. 6 are performed using the same 
parameters, illustrating the reproducibility of the results. 

From the pressures measured at the inlet, middle, and 
outlet of the system, we compute the pressure drops 

AP,{t) = P^nlet{t)-Poutlet{t) and AP„(i) = Prmddle{t)- 

Poutiet{t) (see Figure 2). Figure 4(c) shows a plot of 
these quantities for a typical experiment. After a tran- 
sient characterized by a linear increase, both pressure 
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Figure 4. Temporal variation of imposed and measured expe- 
rimental variables for a typical experiment (Exp. 3, see Table 
I) : (a) Temporal variation of the imposed flow rate (1/15) .Q 
for one syringe. - (b) Ratio of the total injected fluid volume 
to the available pore volume, computed from the imposed flow 
rate using Eq. 2 - (c) Corresponding pressure drop measure- 
ments. Plateaus characterize the different steady-states (ss). 
The sudden drop around t = 75h is an experimental artifact 
due to the reflUing of the syringes, as explained in the text. 
Steady-states obtained before and after refllling are distingui- 
shed using the subscripts a and b. 



drops stabilize and fluctuate around constant average va- 
lues. This behavior is identical to what has been observed 
previously [27] and defines steady state 1 (ssi). The fluc- 
tuations of AP in steady state reflect the dynamics of the 
non-wetting clusters - transport, mergings and snap-offs 
- as the system explores different configurations [27, 28]. 
Immediately after we change the fiow rate to the value 
Q2 > Qi, both pressure drops display a very rapid in- 



crease, and a new steady-state (SS2), characterized by 
higher values of AP, is quickly established. When we 
change the flow rate back to Qi, we observe again a rapid 
variation of the pressure drops towards a third plateau 
defining steady-state 3 (ssa). 

Because experiments are performed at slow flow rates, 
and to obtain enough statistics on the measurements, 
they must typically run for several days. Indeed, the qua- 
lity of the statistics is determined by the ratio Np^, (t) of 
the total volume of fluids injected into the system to the 
total pore volume of the porous matrix, namely : 



Np, (t) 



Q jt) .t 
WLaS' 



(2) 



where Q denotes the total flow rate, t denotes time and 
W , L, a and (j) are the width, length, thickness and poro- 
sity of the porous matrix, respectively. Figure 4(b) shows 
Npy, calculated from the imposed flow rates as a func- 
tion of time, for a typical experiment. In all the experi- 
ments presented here, the duration of the steady states 
correspond to the injection of 0.3 to 1.4 pore volumes 
in the system. As a consequence, it is necessary to re- 
flU the syringes one or several times in the course of an 
experiment. This is systematically performed using the 
following protocol : the outlet and inlets of the model 
are closed and the syringe pump is stopped immediately. 
The reflUing process takes approximately 20 minutes. Be- 
cause the vents used to close the model are not perfectly 
air-proof, the pressure in the system relaxes towards at- 
mospheric pressure during this process, explaining the 
sudden drop of AP observed on Figure 4(c). However, by 
looking at the pictures recorded throughout the process, 
we have checked that this does not affect the flow struc- 
ture, which remains immobile over the duration of the re- 
flUing procedure. Furthermore, we observe that once the 
syringe pump is restarted, AP retrieves its initial value 
after a delay originating from the compressibility of air. 
Thus this refilling procedure does not affect the results 
of the experiments. In the following, when necessary, we 
distinguish steady states obtained before and after refi- 
ling using the subscripts a and b, respectively (see Figure 
4(c)). 



B. History-independence of pressure drops 

Figure 5 shows the temporal evolution of the pressure 
drops APi and AP^ for the 6 experiments listed in Table 
I. In all cases, we observe a behavior similar to what has 
been described in the previous section, namely rapid va- 
riations of AP upon changes of flow rate between pla- 
teaus characterizing steady states. Most importantly, we 
observe that whatever the value and sign of AQ, the pres- 
sure drops are similar in ssi and SS3. This suggests that 
the steady-state pressure drop APgg only depends on the 
imposed flow rate, and not on the history of the system. 



Note also that for a given flow rate, APgg values are re- 
producible from one experiment to the next, regardless 
of how steady state has been reached. 
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Figure 5. Measured pressure drops for the 6 experiments 
listed in Table I. 

To obtain a quantitative indication of this history- 
independence, we compute average steady-state pressure 
drop values (APgs) over the duration of each steady-state, 
thus over periods of time corresponding to the injection 
of 0.3 to 1.4 pore volumes in the system, as already in- 
dicated earlier. For each experiment, we compare the va- 
riation of (APss) between ssi and SS3 to the fluctuations 
of APss within these two states. For this we compute the 
ratio : 



X 



{AP,s,)-{APs 



s3/ 



0.5 {a (APss,) + <J (APss,)) 



(3) 



where (• • •) represents a temporal average over the dura- 
tion of a steady-state and (t(APssi) « ^(APssa) are the 
standard deviations of APgg in ssi and SS3. The values 
of X computed from the inlet and middle pressure drops 
for the different experiments are reported in Table I. As 
can be observed, almost all of these values are smaller 
than one. Slightly higher values are observed in Exp. 3 
and Exp. 4. However, our experimental temperature data 
suggests that temperature-induced viscosity fluctuations 
most likely explain this fact. Indeed, the largest tem- 
perature variations both between ssi and SS3 and wi- 
thin ssi or SS3 occur for these two experiments. There- 



fore, we find that within the precision of our measure- 
ments, the steady-state pressure drop values are history- 
independent. 



C. History-independence of saturation and 
non-viretting cluster size distributions 

Pressure measurements suggest that the observed 
steady-state only depends on the imposed flow rate. Ho- 
wever, these measurements do not give us detailed spa- 
tial information. Thus, we now analyze the images of 
steady-state flow patterns. As explained earlier (see Fi- 
gure 3), the flow structure consists of clusters of the non- 
wetting phase, air, surrounded by the wetting viscous 
water-glycerol solution. As air clusters are transported 
through the porous medium, they are fragmented or mer- 
ged, giving many different realizations of the flow pattern. 
However, within steady-state, the saturation and the size 
distribution of air clusters remain constant on average 
[27]. Thus, similarly to what we have done for pressure 
drops, we compute average steady-state saturations and 
cluster size distributions and compare their properties in 
ssi and SS3. To obtain a good statistics on these quan- 
tities, we have chosen the frame rates so that successive 
images sample different realizations of the flow pattern. 
All averages are computed over series of images, typically 
100, spanning a time range corresponding to the injection 
of 1/5 of the pore volume at the minimum. 

As explained in Section II, the gray-scale histograms of 
the raw images give us a direct indication of the propor- 
tion of the two phases in the system, thus of saturation 
(see Figure 3). Figure 6 shows average steady state histo- 
grams for the 6 experiments listed in Table I. It is clear on 
all these graphs that the histograms corresponding to ssi 
and SS3 can be distinguished from those corresponding to 
SS2 : indeed, as expected, the saturation depends on the 
flow rate, as reflected by different air and liquid "peak" 
heights on the histograms. However, histograms are simi- 
lar in ssi and SS3, suggesting that the steady-state satura- 
tion is also history-independent. It is important to note 
that histograms are directly obtained from raw images 
without prior processing, which excludes eventual arti- 
facts due to image processing. However, they do not give 
us any information about the spatial repartition of the 
two phases in the system. 

To obtain this information, we process the images 
to identify air clusters and compute their size distribu- 
tions. Image processing is performed using ImageJ [39]. 
Raw images are thresholded to obtain binary (black 
and white) images on which we run a standard particle 
analysis algorithm^ to identify air clusters and measure 
their sizes n. From steady-state images, we compute the 
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1. We use the "Particles4" ImageJ plugin written by G. Landini 
(see http ://www.dentistry.bhani. ac.uk/landinig) 



Figure 6. Gray-scale image histograms averaged over steady- 
state images for the 6 experiments listed in Table I. Note that 
no image processing has been applied before obtaining these 
histograms. Data for steady-states ssi , SS2 and SS3 are repre- 
sented in black, blue and red, respectively. Different symbols 
refer to averages performed over different series of 100 images, 
namely : o, ■ for images in ssia, -I-, x in ssu, □, {> in SS2 A, 
V in ss3a and <l, o in ss-it- 



normalized probability density functions of n, i.e. non- 
wetting cluster size distributions {p{n)), where (• • •} re- 
presents an average over a series of w 100 images. 

Figure 7 shows the distributions {p{n)) computed for 
the 6 experiments listed in Table I. These distributions 
typically display a power-law-like behavior with a cutoff 
at large cluster sizes [27, 28]. As mentioned by Tallaks- 
tad et al. [27], the obtained distribution is affected by 
threshold values, which must thus be carefully chosen 
using visual inspection. Here, we focus on the variations 
of the distribution with the history of the system. There- 
fore, the most important requirement is that the image 
processing procedure is used consistently throughout one 
experiment. To avoid possible bias due to variations of 
illumination in the room, the experimental setup is iso- 
lated behind a dark curtain. The camera exposure time 
and aperture are the same for all experiments, and we 
use the same thresholding parameters, carefully chosen 
by visual inspection, for all experiments. This allows us 
to compare images obtained in ssi, SS2 and SS3 for a gi- 
ven experiment and from one experiment to the other in 
a meaningful way. Similarly to what we observed for the 
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Figure 7. Average cluster size distributions {p{n)) computed 
from steady-state images. We use the same symbols and colors 
as on Figure 6 



histograms, it is possible to distinguish the ssi and SS3 
distributions from those corresponding to SS2 (see Figure 
7). This is coherent with the results of previous studies 
indicating that distributions are shifting towards higher 
clusters sizes when increasing the flow rate [27, 28]. Ho- 
wever, the ssi and ssi, distributions are similar, meaning 
that the steady-state non-wetting cluster size distribu- 
tions are history- independent. We have checked that whe- 
reas varying the threshold values affects the distributions, 
typically by shifting then towards lower or higher cluster 
sizes, it does not modify the results in terms of history 
independence. 

The experimental boundary conditions imposed that 
the controlled flow variables were the total flow rate and 
the fractional flows. In the next Section, we turn to nu- 
merical simulations to further investigate the history- 
dependence of the steady-state for different boundary 
conditions, as well as higher Ca values and different vis- 
cosity ratios M . 



IV. THE NETWORK MODEL 

The two-dimensional experimental porous medium is 
modeled by a network of tubes orientated at 45° with 
respect to the overall flow direction. The tubes (or links) 




Overall flow direction 



Figure 8. Illustration of the network formed by tubes that 
are connected to each other through nodes where the dashed 
lines are intersect. One single tube is colored gray. Spherical 
glass beads makes the tubes as hour-glass shaped. 



intersect at the vertices (or nodes) of the network with co- 
ordination number 4. The nodes are considered to have 
no volume, so the tubes consist of both the pore and 
throat volumes. The network is illustrated in Figure 8. 
The disorder due to the random position of the glass 
beads in the experiment is introduced in the model by 
choosing the radius r of each tube randomly from a 
uniform distribution of random numbers in the range 
[0.1Z,0.4/], where I is the length of a tube. In order to 
incorporate the shape of the pores in between spherical 
glass beads, each tube is considered as hour-glass shaped 
which introduces the capillary effect in the system. The 
network transports two immiscible fluids, one of which 
is more wetting than the other with respect to the pore 
walls. The fluids are separated by menisci and we ob- 
tain the capillary pressure pc at a meniscus inside the 
hour-glass shaped tubes from a modified form of Young- 
Laplace law [40, 41], 



27 cos n 
Pc = [1 



27rx 



(4) 



where a; is the position of the meniscus and 7 is the inter- 
facial tension between the fluids. 9 is the wetting angle 
and we consider perfectly wetting conditions, i.e. 6 is 
either 0° or 180°. The flow is driven by maintaining a 
constant total flow rate Q throughout the network which 
introduces a global pressure drop. The local flow rate q in 
a tube with a pressure difference Ap between its two ends 
follows the Washburn equation of capillary flow [40, 42] 



ak 



q = 



Mcff'(Snw)^ 



Ap-^Pc 



(5) 



where k 



is the permeability for cylindrical tubes. 



Any other cross-sectional shape of the tubes will lead 
only to an overall geometrical factor. Here a is the cross- 
sectional area of the tube and /J.off(snw) is the volume 
average of the viscosities of the two phases present inside 
the tube. Hence, it is a function of the local non-wetting 
saturation s^w in that tube. The sum over pc runs over all 
the menisci inside the tube. The property that the fluid 
flux through every node will be zero is used to obtain the 
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local pressures at the nodes. The set consisting of one 
equation (5) per tube, together with the Kirchhoff equa- 
tions balancing the in and out flow at each node are then 
solved using the Cholesky factorization or the conjugate 
gradient method. The system is then integrated in time 
using an explicit Euler scheme. Inside a tube all menisci 
move with a speed determined by q. When a meniscus 
reaches the end of a tube, new menisci are formed in the 
neighboring tubes. In each link, a total maximum num- 
ber of menisci is allowed to form. When the number is 
exceeded, two nearest menisci are merged keeping the vo- 
lume of each fluid conserved. Here, we have considered a 
maximum of 4 menisci in one tube (i.e., 2 non- wetting 
bubbles) , as it is not very likely to form a lot of menisci 
in one pore as seen from the experimental observations. 
Further details of the model and how the menisci are 
moved can be found in [40, 43]. 

In order to reach the steady state in the simulation, 
we considered two different approaches. The conventio- 
nal way is to implement the bi-periodic (BP) boundary 
condition, where the links at the outlet row are connec- 
ted to the inlet links, so that the network acquires a to- 
roidal topology [26]. In this case the simulation can go 
forever, regardless of whether one fluid percolates or not. 
However, in order to keep the flow going, the global pres- 
sure gradient is maintained by considering an invisible 
cut through the system in terms of the pressure. Since 
the system is closed with this boundary condition, the 
individual fluid volumes remain constant throughout the 
simulation. The non-wetting saturation S'nw = Kiw/^ is 
therefore an independent parameter here, along with the 
total flow rate Q, whereas the non wetting fractional flow 
Fnw = Qnw/Q fluctuates over time. 

Implementing the bi-periodic (BP) boundary condition 
is of course impossible in the experiments. As described 
before, in the experimental setup, two fluids are injected 
at one edge of the Hele-Shaw cell through a series of al- 
ternate inlets and the opposite edge is kept open. In this 
case, the flow rates of the two fluids can be controlled 
independently. Thus, the control parameters are the to- 
tal flow rate Q and the fractional flow Fnw, whereas the 
saturation S^w fluctuates. In order to have a close emula- 
tion of the experimental ensemble, we also implement the 
open boundary conditions (OB) in the simulations, where 
the individual flow rates at the inlet links are controlled. 
Therefore, in OB, the system is open in the direction of 
total flow while we consider periodic boundary conditions 
in the direction perpendicular to the overall flow. 



V. SIMULATION RESULTS 

Simulations are performed considering networks of 
256 X 256 hnks for BP and 128 x 192 hnks for OB. In 
order to avoid any traces from the inlets in OB, only a 
128 X 128 segment of the network towards the outlets is 
considered for the analysis (see Fig. 10). Each link has 
a length of 1 mm, which is equal to the bead diame- 



ter used in the experiments. Most of the simulations are 
performed for the viscosity ratio M = 1. Three different 
capillary numbers, Cai — 1.92 x 10~^, Ca2 = 9.15 x 10~^ 
and Gas = 2.88 x 10"^ are considered for BP. For OB, 
the capillary numbers considered are Cai = 3.2 x 10^"^, 
Ca2 = 3.2 X 10-2 and Gag == 7.0 x lO-^. For OB we 
choose a similar fractional flow F,^„ = 0.5 as in the expe- 
riments. In BP, we run the simulation for the saturation 
Snw = 0.74 which we find close to the critical point for 
the range of parameters we considered here [32] . We will 
report only one set of simulations for M = 10""' with 
BP for a system with 128 x 128 links and Snw = 0.5 
with capillary numbers 5.06 x 10"^ and 1.24 x 10"^, as 
the conjugate gradient solver converges very slowly for 
viscosity unmatched fluids, making the simulation very 
computationally expensive. 

To investigate the history dependence, we use a proce- 
dure analogous to the experimental one (see Sec. Ill A) : 
the simulation is started from the initial condition with 
a capillary number setting a constant flow rate Qi. Once 
a steady state ssi is reached, the capillary number is al- 
tered to a different value setting another constant flow 
rate Q2 , and this new flow rate is maintained until a dif- 
ferent steady state SS2 is reached. Finally, the capillary 
number is set back to the initial value and the system is 
allowed to evolve towards a third steady state SS3. Each 
simulation thus involves two different capillary numbers 
among Gai, Ga2, and Gag, therefore 6 different simula- 
tions have been performed for each boundary condition. 
The steady states are compared with three different cri- 
teria - the average pressure drop, the distribution of fluid 
saturation over the system, and the non-wetting cluster 
size distribution. All the measurements are averaged over 
50 to 100 configurations in the steady-state and 5 to 10 
different realizations of the network. 

Fluid morphologies for a typical simulation in the three 
steady states ssi, ss2 and SS3 are shown on Figures 9 and 
10 for BP and OB respectively. Here the simulation starts 
from Ca = 9.15 x 10^^ to reach ssi, then it is altered to 
Ca = 2.88 X 10^2 to reach SS2, and then it is again turned 
back to the initial Ga = 9.15 x 10"^ to reach SS3. Figures 
9 and 10 show the distribution of saturation over the 
network in these three steady-states in (a), (b) and (c) 
respectively, where the gray-scales from black to white 
correspond to Snw = 1 to inside a link. For BP, it is 
not possible to see any difference in the gray-scale satu- 
ration distributions, as the system is closed and the total 
saturation is conserved. However, as we will see, iden- 
tifying the clusters allows us to distinguish ssi and SS3 
from SS2 (sec Figs. 9 (d), (e) and (f)). In OB, the satu- 
ration distributions look similar in ssi and SS3 whereas 
SS2 shows more non-wetting saturation than the others. 
This is consistent with previous numerical studies sho- 
wing that the variation of saturation with Ga depends 
on the viscosity ratio, fractional flow and other flow pa- 
rameters [26]. We then identify the non- wetting clusters 
using Hoshcn-Kopclman algorithm [44] . As every link can 
be occupied with both the fluids, a clip threshold in the 
















Ca = 9.15 X 10-^ Ca = 2.88x10"^ Ca = 9.15x10"^ 

Figure 9. Typical steady-state fluid morphology over the net- 
work in BP for M = 1 and Snw = 0.74. The distribution of 
non-wetting fluid saturation inside the links in the steady- 
states ssi, SS2 and sss are illustrated in (a), (b) and (c) res- 
pectively and the corresponding capillary numbers are indica- 
ted under each column. The gray-scales from black to white 
correspond to the non-wetting saturation from 1 to inside 
a link. The non-wetting clusters identified by the Hoshen- 
Kopelman algorithm are shown by different random colors in 
(d), (e) and (f). 



link saturation is considered to identify the clusters [45] . 
If a neighboring link has a non- wetting saturation higher 
than the clip threshold, the link is then considered to be- 
long to the same cluster. The clusters are shown in (d), 
(e) and (f ) of Figures 9 and 10 for the three steady-states, 
where each cluster is drawn in a different color, chosen 
randomly. The distribution of the clusters shows a clear 
characteristic difference of ssi and SS3 from SS2 both in 
BP and OB. The cluster sizes in ssi and SS3 look very 
similar, and they are distinctly different from that of SS2. 

The temporal evolution of pressure drops for different 
simulations is shown in Figures 11 (M = 1) and 12 
(M = 10^4) for BP and in Figure 13 for OB. The two 
different capillary numbers for each simulation correspon- 
ding to ssi, SS2 and SS3 are indicated in the plots. In BP, 
we measure the global pressure drop AP over the whole 
system. In OB, we measure the pressure drops Ai-^ at 
the inlet nodes and AP^, at the middle of the system, 
with respect to the outlet where the pressures are ave- 
raged over all the nodes in the corresponding row, in 
the direction perpendicular to the flow. We observe very 
similar behaviors in the pressure curves as in the expe- 
riments. In the case of BP, the system is initialized by 
filling the tubes with the fluids randomly at the desired 
saturation Snw, which will be constant throughout the 
simulation. This random initialization has the advantage 
of decreasing the simulation time significantly, since the 
steady state is reached faster than with an initial condi- 
tion in which the two fluids are completely segregated. 
Due to this initial random filling in BP, the global pres- 
sure AP starts from a higher value and decreases with 
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Ca = 9.15x10"^ Ca = 2.88x10'^ Ca = 9.15x10"^ 

Figure 10. Steady-state fluid morphology over the network 
in OB for M = 1. The flgures are drawn in the similar way as 
that of BP. 



time due to the formation of clusters. Subsequently, it 
reaches the steady state ssi and AP fluctuates around a 
constant average value as seen on Figures 11 and 12. In 
OB, the system is initialized by saturating the network 
completely with the wetting fluid and then the simulation 
is started by injecting two fluids simultaneously through 
a series of alternate inlets. The flow rates of individual 
fluids are controlled to obtain the required fractional flow 
Pnw Both drainage and imbibition therefore take place 
at the pore level creating new menisci, which increase the 
pressure drop with time as seen in Figure 13. Away from 
the inlets, the trace of the injection channels vanishes, 
and a steady-state ssi is attained, and APj and AP^ 
fluctuate around constant average values. Next, as soon 
as the capillary number is changed to a different value, a 
rapid change in the pressure drops is observed for both 
BP and OB, and the new steady-state SS2 is reached, 
characterized by different constant values in the average 
pressure drops. When the capillary number is again alte- 
red to the initial value to reach the steady-state SS3, we 
find that the global pressure drops change back to the 
initial average value. Moreover, the global pressure drops 
corresponding to the same capillary numbers in different 
simulations have the same average value, no matter from 
which condition they have been reached. 

The global pressure estimates therefore completely 
support the experimental observations, i.e. that the 
steady state only depends on the imposed flow rate, and 
not on the initial condition. Now, in order to find a de- 
tailed microscopic information in this regard, we mea- 
sure the distribution of link saturation over the system. 
This measurement provides us with similar information 
as that of experimental gray-scale image histograms, des- 
pite being computed slighly differently. More precisely, in 
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Figure 11. Global pressure drop AP as a function of time 
for BP in different simulations with M = 1 and Snw ~ 0.74. 
A rapid change in AP can be observed as soon as the overall 
flow rate is altered. When the flow rate is restored to the 
initial value, AP settle back to the initial steady-state value 
as seen for ssi and SS3. 



Figure 13. Time evolution of global pressure drops APi at 
the inlet and APm at the middle of the system for OB with 
M = 1. A rapid change in both the pressure drops can be 
observed as soon as Ca is altered, but they again settle back 
to the initial value when the flow rate is restored to the initial 



the experiment, the gray-scale of each pixel is counted, 
where one pixel corresponds to any of the three compo- 
nents - the viscous fluid, air or the glass beads. In the 
simulation, on the other hand, we count the non-wetting 
saturation inside each link and compute the histogram of 
the link counts. Therefore, one should not try to make 
a direct match of the histogram patterns from the expe- 
riments to the simulations. The histograms in the three 
steady-states ssi, SS2 and SS3 in different simulations are 
plotted on Figures 14 (M = 1) and 15 (M = lO"") for BP 
and 16 for OB. In each simulation, it is very clear that 
the histograms for ssi and SS3 fall on each other whereas 
they are distinctly different from that of SS2. Moreover, 
the histogram patterns corresponding to the same Ca 
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Figure 12. Global pressure drop AP with time for BP with 
M = 10-* and Snw = 0.5. 



in different simulations are identical. A minor difference 
in the histograms for ssi and SS3 is observed only for 
Ca = 1.92 X IQ-^ in BP. This is due to the appearance, 
at low Ca and high saturation in BP, of percolating non- 
wetting flow channels yielding a high non-wetting frac- 
tional flow {Fnv, ~ 0.99), as shown in the inset of Figure 
14 (e), while the rest of the system is immobilized. 

On the histograms, we observe two distinct peaks in 
BP for Ca = 9.15 x 10^^ : one at Snw ^ 0.8 corres- 
ponding to the links mostly filled with non- wetting fluid, 
and the other at Snw ^ 0.4 corresponding to the links 
mostly filled with wetting fluid. Therefore links can be 
divided into two categories : either highly saturated with 
the non-wetting fluid or with the wetting fluid, rather 
than containing a mixing of two phases. This in turn indi- 
cates the presence of large clusters at this Ca, as already 
observed in the fluid morphology on Figures 9 (d) and 
(f). For the other two Ca values in BP, the histograms 
arc changing towards having one peak, with a more flat 
shape. For M = 10^"* in BP, and in OB, the histograms 
also display one peak, roughly centered but whose posi- 
tion shifts along the x-axis with Ca. This indicates that 
most of the links are filled with a mixture of both fiuids. 
Therefore it is difficult to obtain large clusters in such 
conditions, as observed in the fluid morphology for OB 
(Figure 10). 

Finally we compute the non- wetting cluster size distri- 



11 




Link saturation 



Linlc saturation 



— sS|:Ca = 9.15xlO" 

- - sSj! Ca = 2.88x10"' 
■-■ ss ■Ca = 9.15xlO"' 



(c) 




— sS|:Ca = 2.88x10" 
-- sSj:Ca = 9.15x10"^ 
■-■ ss,:Ca = 2.88x10"' 



(d) 



Link saturation 



Link saturation 

















r^ 


f 


— sS|:Ca= 1.92x10 ' I 


" 0! 






-- ssj: Ca = 2.88x10"' | 


I 1 


■-■ ss,:Ca= 1.92x10"' | 




1,1 




. 




" '' ^^ [ 


/"^v /'' "* '■ 





— sS|:Ca = 2.88x10"' 


(f) 


_ -- ss,:Ca= 1.92x10"' 




■-■ ss,:Ca= 2.88x10" 


/ \ 


- 


/ \ 


,— , 


/ (- 


/' ""'^-^^ y 


1 









(a) 


— SS| :Ca- 3.2x10" . 


- 


'^ -- ss^:Ca = 3.2x10"' - 






\ - ss^:Ca = 3.2x10"^ - 


- 


//' 


'"" \^^N : 


r 


// 


>^ 


'w 







Link saturation 



(c) 


— sSj :Ca= 3.2x10"" . 


-- ss,:Ca = 7.0x10" - 




■-■ ss,:Ca = 3.2x10" - 


- 


/''-^ 


^ 


^T'^Vn 


/' 


\>N 


/,- 


%\ 







Link saturation 



(b) 



SS| :Ca= 3.2x10" 
ss,:Ca = 3.2x10""' 
ss, :Ca = 3.2x10"" 




Link saturation 



(d) 



SS| :Ca = 7.0x10" 
ss,:Ca= 3.2x10"" 
ss, :Ca = 7.0x10"" 



Link saturation 





Link saturation 



Link saturation 



Link saturation 



Link saturation 



Figure 14. Histograms of the network links according to non- 
wetting fluid saturation inside the links over the steady-state 
configurations for BP with M = 1 and S'nw = 0.74. In all the 
simulations, histograms are found similar for ssi and sss. Two 
distinct peaks are observed for Ca = 9.15 x 10~^, implying 
that the links are either highly saturated with the non-wetting 
fluid or the wetting fluid. The inset of (e) shows non-wetting 
fractional flow i^nw for the corresponding simulation where 
Kiw ~ 0.99 at Ca = 1.92~^, which is reason for minor varia- 
tion in the histogram patterns for ssi and sss. 



buttons at different steady states. Here, the size n of a 
cluster is defined by the total number of links that be- 
long to the cluster. The probability p(n) to have a n-sized 
cluster is then defined as p(n) ~ N{n)/Ntot, where N{n) 
is the number of n-sized clusters out of total Nfot clusters 
identified. p{'n) is averaged over different configurations in 
the steady state and different samples of the network. In 
Figures 17 and 19, {p{n)) is plotted in log-log scale for BP 
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Figure 15. Steady-state saturation histograms for BP with 
M = 10"* and Snw = 0.5. 



Figure 16. Histograms of the links according to the non- 
wetting fluid saturation inside the links over the steady-state 
configurations for OB with M = 1. In all the simulations, 
histograms are found identical for ssi and ssa. Histograms 
contain only one peak around the middle, which means that 
the links are mostly saturated with the mixture of both the 
fluids. 



and OB respectively, for different simulations. For all the 
simulations, we find that the cluster size distributions are 
identical for ssi and SS3 whereas they are different from 
SS2. Distributions for the same Ca for different simula- 
tions are also the same, showing that the steady-state 
cluster size distributions are history independent. 



VI. CONCLUSIONS 

In this article we have considered the question of 
history dependence in the steady-state two-phase flow 
in porous media and presented detailed experimental 
and numerical investigations in this context. Experimen- 
tally, a quasi two-dimensional laboratory model consis- 
ting a Hcle-Shaw cell filled with glass beads is conside- 
red, through which two phases, a gas-liquid pair with a 
viscosity ratio 10"'* flows at constant flow rate. The sys- 
tem is allowed to evolve to a steady-state where the global 
pressure drop fluctuates around a constant average value. 
Steady-states corresponding to the same control parame- 
ters (e.g. capillary number) are attained from different 
initial conditions. In order to characterize the complex 
flow patterns in the steady-state, the gray-scale histo- 
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Figure 17. Steady-state non-wetting cluster size distribu- 
tions, (p(n)) vs. n for BP with M = 1 and Snw = 0.74. For 
all the simulations, ssi and sss are found to have the same 
distributions, which depend only on the capillary number. 



Figure 19. Steady-state non-wetting cluster size distribu- 
tions, (p(n)) vs. n for OB with M = 1. Similar to BP, the 
distributions are found identical in ssi and ssa in all the si- 
mulations. 



gram of snapshots and the non- wetting cluster size distri- 
butions are analyzed. The system is then modeled nume- 
rically by a network of disordered pores transporting two 
immiscible fluids. Steady-state situation in the network 
model is attained implementing two different boundary 
conditions, the toroidal one with constant saturation and 
the open boundary with constant fractional flow, similar 
to the experiments. The global pressure drop, distribu- 
tion of non-wetting pore saturation over the network and 
the cluster size distributions are computed. Both the ex- 
perimental and numerical results show that when both 
the fluids are flowing in the steady state, different mea- 
surements corresponding to the same control parameters 
are identical, no matter how the steady-state has been 
reached. Thus, unlike the transients, the steady-states 



only depend on the external parameters, but do not de- 
pend on the initial preparation of the system or the his- 
tory of the process. We therefore conclude that, within 
the range of parameters explored in this study, the steady 
states in simultaneous flow of two phases through porous 
medium are history independent. 
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Figure 18. Steady-state non-wetting cluster size distribu- 
tions, (p(n)) vs. n for BP with M = 10~* and Snw = 0.5. 
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